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We describe inhomogeneities in a ACDM universe with a gradient series expansion and show that 
it describes the gravitational evolution far into the non-linear regime and beyond the capacity of 
standard perturbation theory at any order. We compare the gradient expansion with exact inho- 
mogeneous ALTB solutions (Lemaitre-Tolman-Bondi metric with the inclusion of a cosmological 
constant) describing growing structure in a ACDM universe and find that the expansion approxi- 
mates the exact solution well, following the collapse of an over-density all the way into a singularity. 

I. INTRODUCTION 



The study of cosmological inhomogeneities through perturbation theory has been one of the cornerstones of modern 
cosmology [1]. Indeed, the confrontation of theory with increasingly precise observations of the CMB and the large 
scale clustering of galaxies through the extended use of linear perturbation theory has played an important role in 
the emergence of the concordance ACDM cosmological model and has allowed meaningful speculation about the very 
early universe [2]. At higher order, perturbation theory is also being used to probe subdominant non-linear corrections 
and the associated non-Gaussianity, promising to shed additional light on the primordial universe. 

\.J Although crucially important for following the development of cosmic structures in the primordial universe, and 

even more recent eras on large enough scales, standard perturbation theory is naturally limited. In the matter era 
perturbations grow, the density contrast becomes non-linear and the inhomogeneities separate themselves completely 
from the background expansion; at this point standard perturbation theory breaks down completely at all orders. Some 
recent reformulations can somewhat extend its range of validity in the quasi-linear regime [3] but the most complete 
quantitative description of gravitational clustering is necessarily achieved by Newtonian N-body simulations. However, 

-4— > analytic approximations are important for obtaining insights into the process of structure formation. In particular, 

the Zel'dovich approximation along with a range of improvements is a prominent analytic approximation scheme and 

i— —i can be used to understand various features of the resulting cosmic web [4] . 

In this paper we use an expansion different from well known and developed cosmological perturbation theory. We 

►C| approximate the metric perturbations of a ACDM universe by a series whose terms contain an increasing number of 

spatial gradients of the initial metric and not powers of the metric functions (or the density), with each term multiplied 

^^ by a specified function of time. At leading ("zeroth") order this is just the "separate universe" approximation. 

£■ — , However, the inclusion of higher order terms allows for the subhorizon evolution to be traced as well and already at 
the first non-trivial order the approximation goes beyond standard perturbation theory and can describe the non- 

~S^. linear regime with gravitational collapse/void rarefaction. The inclusion of even higher orders approximates the true 

evolution to an increasing degree, at least if followed for a characteristic timescale. It is interesting to note that the 

gradient expansion as will be used below can reproduce the Zel'dovich approximation for the density contrast but 

also allows for relativistic corrections to be included [5, 6]; the latter are indeed crucial for increased accuracy during 

£>. the later stages of gravitational collapse. 

The gradient expansion can be used to evolve any initial configuration of perturbations, for example those described 
by a post-inflationary scale invariant spectrum. In this work, after presenting the general formulation, we compare 
the first two non-trivial terms in the expansion with exact Lemaitre-Tolman-Bondi [7-9] solutions with the inclusion 
of a cosmological constant (ALTB) from Ref. [10] in order to assess the accuracy of the gradient series approximation. 
We find the approximation to provide a fairly good description, certainly much beyond the capacity of standard 
perturbation theory. We include terms with up to 4 spatial gradients and the inclusion of higher order terms is 
expected to fare even better in accurately describing the non-linear gravitational evolution. 

II. THE INHOMOGENEOUS ACDM UNIVERSE AS A SERIES IN SPATIAL GRADIENTS 

The gradient expansion as a method for obtaining solutions to the Einstein Equations dates back to [11] and has 
been used as a means for approaching the initial singularity of an inhomogeneous universe [12]. It was discussed in 
more general terms in [13, 14] and more recently applied to inflation in [15, 16]. A leading order covariant computation 
can be found in [17]. It as also proved useful in studying the backreaction of cosmic structure formation, see [18] and 



more recently [19]. In most of past works the gradient expansion has been applied directly on the Einstein equations. 
However, in [5, 6] a Hamilton- Jacobi (HJ) approach was used which arguably offers a less laborious way to perform 
the expansion. In this section we follow the HJ methods to obtain a gradient expansion for inhomogeneities in a 
ACDM universe. 

The applicability of a gradient expansion does not rely on some quantity being a small perturbation around a 
given background but would generally require the existence of a hierarchy between spatial and temporal derivatives 
d t Q > diQ ~ -J£ > —H where Q is either a matter or a metric quantity with At and L the characteristic time-scale 
and length-scale of the system respectively. One would generically expect the gradient expansion to be applicable 
on length-scales L ^> At i.e. scales beyond causal contact. This has indeed been used in early universe cosmology 
where dtQ ~ HQ under the guise of the super- horizon (-W < x) "separate universe" approximation which is nothing 
but the lowest order in a gradient expansion [20] . However, the gradient expansion can also be used for length scales 
within causal contact provided that the spatial variation of the initial data is smooth enough and the spatial gradients 
small. One could then still use the gradient expansion but for a limited time span: bigger spatial variations would 
correspond to shorter times for which the gradient series would apparently converge. We will make these statements 
more precise below. 

A. Hamilton-Jacobi for ACDM 

Let us start by writing an action for gravity with a cosmologies! constant A and Cold Dark Matter (CDM) 



/ d4xV ^ L ( Wi? " 2A ) " \ p fo^x^x + l) 



S = / d 4 xV=i 



(1) 



where x is a potential for the 4-velocity of CDM, U^ = —g^d^x, an d p is the energy density. The latter acts as 
a Lagrange multiplier whose variation ensures that U^U^ = — 1. Variation with respect to x gives the continuity 
equation Vu(p9 M x) = while variation of the CDM part of the action with respect to g^ v gives the usual CDM energy 
momentum tensor T M „ = pU^U v . Gravity is described by the standard Einstein-Hilbert term with the addition of 
a cosmological constant with ^R the 4-dimensional Ricci Scalar, and K = 8ttG. The ADM decomposition for the 
metric can now be used to develop a Hamiltonian formalism for this system. Writing the metric as 

.goo = -A^ 2 + hijrfN* , g 0i = 7 tf JV J ' , g tj = ^ , (2) 

and defining the canonical momenta 
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T=py/ly]l + r'd i xd i x, (4) 



with 



we can write the action as 



Etj = ~7« - V (< ify , E = h ij E i:j , (5) 
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/ / - ^ (x^ - y ) - ^ (R - 2A) + tt* yJi+jiJdixdjx , (7) 



Hi = -2V fc < + TT x d iX ■ (8) 



Choosing the velocity potential x to define the time hypersurfaces of the system implies that diX — along with 
fit = 1 and N = 1. The metric then takes the form 

ds 2 = -dt 2 + -fij (t, x)dx i dx j , (9) 



and the spatial coordinate lines comove with the matter - there is no peculiar velocity. This choice then corresponds 
to the comoving synchronous gauge 1 . It also means that ir x now plays the role of the Hamiltonian density:— ir x = %. 
By imposing the energy constraint W = 0we obtain in this gauge a hamiltonian density for the gravitational field 

2 « _ ki ( . . . i. . \ \py, 



H = — t^tt*' n k ^i - -7^*1 )-^-(R-2A), (10) 

where R is the Ricci scalar of the 3- dimensional time hypersurfaces. The action finally takes the form 

S= J dtd 3 x L ij ^-U + 2N& k iP\ , (11) 

and defines a constrained Hamiltonian system where the canonical momentum 7r y is constrained to be covariantly 
conserved 

V k Tr ki = . (12) 

Let us now apply the HJ approach to the Hamiltonian system (11) - see eg [22]. The action can be considered as a 
function of the metric and time S[t, 7^]. Then the momentum is simply the variation of this action with respect to 
the metric it 13 = SS/Sjij , the Hamiltonian is minus the rate of change of the action with time and the HJ equation 

dS f ,, f 2k SS 5S ( 1 \ Ji tri „ A .l „ 

contains all the dynamical information of the system. It is a single partial differential equation for S as a functional 
of 7^ and a function of t. Once S[t, 7^] is determined the metric can be obtained from 

-^f = —^s- — Qlikljl ~ lijlkl) ■ (14) 

dt ^7 d-fki 

Furthermore, the remaining constraint (12) will be automatically satisfied if 

S= I ' d z x^F{t, lv] ), (15) 



where T{t, 7jj) is some scalar function of the metric making S invariant under 3-D diffcomorphisms. Indeed, the 
variation of such a functional with respect to the metric will yield a covariantly conserved tensor which will thus 
satisfy (12). Once a solution is obtained, the CDM density can be obtained from 

/3 (t,x) = — 4° (16) 

where po is a constant. Equations (13) and (14) are the basic equations that need to be solved in this approach. 

B. The gradient expansion 

The gradient expansion in this formulation begins by writing an ansatz for J 7 as [5, 6] 

T = -2H(t) + J{t)R + Li(t)i? 2 + L 2 (t)R lJ R l3 + ... , (17) 

a series in powers of the 3-D Ricci curvature (involving an increasing number in gradients of 7y). The HJ equation 
(13) can now be solved separately for terms with a different number of gradients. This gives for the time dependent 



1 There is conceptual similarity of this choice to the Lagrangian picture of cosmological perturbations [21]; see also the discussion below 
on the relation to the Zel'dovich approximation. 



coefficients in (17) 
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e.t.c. The solution for H then determines all other functions. For A ^ equation (18) gives 



/A , /V3A x 



(22) 



As the notation anticipated, H is simply the background Hubble parameter. The other functions arc then given by 

(sinh u) du , 
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where 



L(t) = e 1 *" 
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ti is the initial time and u — ^-|^ t. Given the above, the metric can now be obtained by solving equation (14). Up 
to terms with 4 gradients we obtain 

^ = 2H 7ij + J{R lij - AR tj ) 

-\-L\ (3^/ijR — SRRij + 8i?ijj) 
+L 2 (z lij R kl R kl - m ik R k 3 + 3 7lJ i?l fc | fe 
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So far, the resulting Eq. (27) is exact if one includes all orders. The gradient approximation now consists of solving 
this equation iteratively for 7^ . At zeroth order in gradients we have 
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which gives 



1 ( ff{t) = A\t)k ij , 



(29) 



where kij is the initial metric at time t = ti and 



A(t) = e 



-X h - 
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This is simply the separate universe approximation where each point evolves as a homogeneous FRW universe. Up 
to second order in gradients we get 



5 7 (2) 
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where we have substituted the order result in the terms that contain two gradients: Rij = Rij(ki m ) and R = R(ki m ) 
are the Ricci tensor and Ricci scalar of the initial metric kij. The solution of (31) gives 



jV = A 2 (t)k l0 + XQkXRkij - 4R tJ ) 



*3 ) i 



where 



(32) 



\(t) = e 2f t" 



~ 2 ft H T 

e Jt i J 



dx = 



V3A 



. , V3A 

Sinn 1 

2 



V3A + 
4/3 2 l 



(sinhu) ' J{u)du. 



(33) 



Proceeding to 4 gradients we have 

a 7 (4) 
at 
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-D\R) ikhj + D 2 R\ij + D 3 Rij' \ k , 



(34) 



where a vertical bar denotes covariant derivation with respect to the seed metric kij and the Bianchi identity has been 
used to eliminate some terms. The C and D coefficients then read 
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We obtain 
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{R ifcfejj — 4Rij> \k + R\ij 
One could proceed in this manner to obtain higher order terms. 2 



(36) 



C. Range of validity and the relation to the Zel'dovich approximation 



Before proceeding further it is worth estimating the range over which the gradient expansion can be applied. 
Focusing at early enough times the cosmological constant doesn't play a role and the universe evolves as if it were 
matter dominated. The metric in the gradient expansion then reads 
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2 A gradient expansion solution equivalent to (36) was also obtained in [23] using a method different from the Hamilton-Jacobi approach 
utilized here. It would be interesting to further explore the relation with that approach. Our thanks to the referee for pointing it out. 



where we have ignored terms that are sub dominant when t ^> t\. The 2-gradient and 4-gradient terms become of the 
same order of magnitude as the Oth order one after a time t 



O(few) 



1 



t 2 R?l 2 



or t c 



O(few) 



t 2 (V 2 R) 3 / 4 



(38) 



The exact coefficient is of course dependent on the form of the initial 3- metric kij . We can therefore estimate that the 
expansion should be valid for t < t con . After this time the contributions to the metric of successive terms containing 
more gradients are no longer perturbatively ordered and the iterative solution to (27) breaks down. However, note 
that the timescale over which the series breaks down is commensurate with the time of collapse for a region of positive 
spatial curvature (where Rkij — ^Rij < 0). After that the region either forms a singularity or a caustic is created 
which would need extra physics for its resolution. Either way, one would not expect the current description to hold 
anymore and thus the limited timescale over which the expansion holds does not appear to be a serious limitation. 
For negative curvature regions which keep on expanding things are less clear-cut since there is no natural limit on the 
evolution timescale. Nevertheless, it would seem that the gradient expansion can provide a fair description for times 
longer than t con even in this case - see [19]. 

Let us now apply the above to our universe. Assuming the standard inflationary initial conditions, the initial seed 
metric takes the form 
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where $(x) is the primordial Newtonian potential and we have scaled the seed metric such that in the absence of 
perturbations the scale factor would be unity today. Note that the lowest order in the above expansion is simply the 
separate universe approximation where each spatial point scales as a FRW, CDM, universe with scale factor a ex t 2 / 3, . 
With kij as the initial metric, we have 
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where 
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The timescale for which the above approximation for the metric is accurate is estimated to be 
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in accordance to (38). Equivalently, applying the Poisson equation V 2 $ = ^a 2 H 2 Q m 5 m , we find an upper bound for 
the matter overdensity S m = Sp m /p m today 



JmO 



1.5/n* 



(43) 



We see that for describing an inhomogeneous universe which resembles ours for the whole of to, the metric can include 
perturbations down to about k = 0.3/iMpc _1 , or 5 m o — 5. Of course, even shorter scales can be accurately described 
but for shorter times, comparable to the collapse times of the over-dense regions. 
Given these approximations, the local density of matter can be obtained from 
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Note that if this expression is expanded to linear order in $ one recovers the result of linear perturbation theory for 
the density contrast in the synchronous gauge. Dropping the B^ term we have 
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which reproduces the well known Zel'dovich approximation. Thus, the gradient expansion in the comoving synchronous 
gauge can be thought of as providing a relativistic extension of the Zel'dovich approximation [5, 6] 3 . In fact, as will 
also be demonstrated below, the higher order terms are crucial for increased accuracy during the more advanced 
stages of the gravitational evolution. 

At late times H — > J ^ and the universe enters a de Sitter phase in which case it is easy to see that the 4- 

gradient terms in the metric are exponentially suppressed, the 2-gradient terms freeze while the O-gradient terms are 
exponentially enhanced. Thus, the gradient expansion becomes increasingly better with time, converging with time 
to the separate universe picture as would be expected, apart from patches which suffer gravitational collapse before 
the de Sitter phase. 

III. COMPARISON WITH ALTB 

A. The LTB metric 

Now let us assess the quality of the Gradient Expansion by comparing it to an exact metric that approximates the 
formation of structure in a ACDM universe. The Lemaitre-Tolman-Bondi metric is spherically symmetric, and can 
describe a central density perturbation embedded in an exact Friedman-Lemaitre-Robertson- Walker (FLRW) metric 
by choosing an appropriate density profile. With the inclusion of dust, curvature and a cosmological constant this 
hence describes the rarefaction of a void as well as the initial expansion and later collapse of an over-density in a 
ACDM universe 4 . The metric is given by 

ds 2 = dt 2 -X 2 (r, t)dr 2 - Y 2 (r, t)dfl 2 , (46) 



with X{r, t) = d r Y(r, t)/ \/l + 2E(r). The metric is described by three free radial functions, being the curvature E(r), 
the matter distribution p(r) and the Big Bang time £bb(?~)- -Due to the gauge freedom in the radial coordinate r, one 
of the three functions can be set arbitrarily, and we choose a gauge in which L Y'(r,t)Y 2 (r,t)p(r,t)dr oc r 3 . Therefore 
the choice of curvature profile E(r) and bang-time profile £bb(t) defines the metric completely. Moreover, now any 
deviation of £bb(?") from a constant corresponds to a purely decaying mode. Following the inflationary paradigm we 
know that the decaying mode must be small such that is has decayed today, and we can safely ignore it by setting 
^BB( r ) = everywhere. All freedom now hence lies in E(r). 
The Einstein equations simplify to 
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(47) 



where Y(r, t) = ra(r, t), and we see that for E(r) oc r 2 we obtain exactly the FLRW solution. Moreover, if we choose 
E(r) such that beyond a certain radius, say r = L, E(r) oc r 2 , while below that radius it has a different shape, that 
E(r) describes an inhomogeneity embedded in FLRW. We choose, 

E(r) = -r 2 kW 7 (^,0] , (48) 



L 

where k = ± 87r ° and the function W n (r/L, 0) is defined in Ref. [10], and interpolates from W n (0, 0) = 1 to W n (l, 0) = 
while remaining C n everywhere, such that the metric is C n . By choosing n — 7 we are guaranteed to have all 
functions entering the problem to be at least C° everywhere. Having W n (0, 0) = 1 and W n (l, 0) — 0, we have defined 
a curvature profile with non-zero curvature at the center of the metric at r = 0, and zero curvature at finite radius 
r = L such that the metric there becomes the exact spatially flat FLRW metric with dust and a cosmological constant 
for r > L. The magnitude of k is arbitrary, and we chose the value it has such that it describes a central over-density 
that reaches the singularity today, and an under-density that is non-linear today but does not suffer form shell crossing 
yet. 

We use the full solution to the Einstein equation as well as Y'(r, t) from Ref. [10] valid for any distribution of dust 
and curvature, until shell crossing occurs. In calculating the higher order gradient terms we encounter second spatial 
derivatives of the spatial Ricci scalar, which contain both Y"(r,t) and Y'"(r,t). Since we have exact solutions to 
Y(r,t) and Y'(r,t), we can safely evaluate the second and third radial derivative of Y(r,t) numerically. 



3 Planar solutions in the presence of A akin to Zel'dovich pancakes are described in [24] and [25] for GR and Newtonian Gravity respectively 

4 Generalizations to spheroidal configurations can be found in [26]. 
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FIG. 1: The radial component of the metric y/—g rr = X(r,t) as a function of radius rata number of times, where different 
colours correspond to different times. The solid lines are the exact solutions, the coloured dots are the gradient-expansion 
solutions, up to fourth order in gradients. On the left we show the solution for an over-density, compensated by an under-dense 
shell and matching to FLRW at exactly r = L. The center of this particular configuration reaches the singularity at t = f 3.5 Gyr 
after the initial singularity. On the right we show an under-density with the same curvature as the over-density, but opposite 
sign. By eye, both scenarios are very well approximated by the gradient expansion. See Figures 2 and 3 for quantification of 
this statement. 
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FIG. 2: The radial component of the metric y/—g rr — X(r,t) at r — 10 L (practically at the centre) as a function of 
time t, comparing the exact solution (solid red), linear perturbation theory (dotted black, uppermost line), 2nd order gradient 
approximation (dotted green, innermost curve) and 4th order gradient approximation (dotted blue). While linear perturbation 
theory clearly does not show any form of collapse, both lowest (2nd) and next order (4th) in the gradient approximation show 
the collapse. The inclusion of 4th order gradient terms approximates the exact solution better than the truncation to 2nd 
order. 



B. Comparison of the gradient expansion with the exact solution 



We calculate the initial conditions to Eq. (27) from the exact metric satisfying (47) at some early initial time and 
numerically integrate the set of six coupled differential equations that describe Ait), J(i), L(t), X(t), J X(t) J (t) / A 2 (t) dt 
and J L(t)/A 2 (t) dt at any later time. We then compare the resulting solution (36) to the exact solution from (47). 

In Figure 1 we show the results for an over-density with k — „ ° (left) and an under-density with k = g— - 

(right). We plot y/—g rr = X(r, t) as a function of radius r for different times, where time is encoded by the colouring 
of the lines and dots. Solid lines correspond to the exact solution, while dots present the expansion up to fourth order 
in gradients. The expansion follows the features in the metric very well, in both the compensating shells and the 
central over- and under-densities. Obviously, for r > L, where the metric is exactly FLRW, the gradient expansion 
corresponds to the exact solution with g rr ex a 2 (t). 

To clearly show how well the gradient expansion follows the gravitational collapse of an object, we show in Figure 2 
the radial metric component \J—g rr = X(r,t) very close to the centre. Exactly at the centre we have X(r,t) = OVi, 
since it corresponds to r a(t) in FLRW, which is why we plot X(r,t) just off-centre, at r = 10 _5 L. For illustration 
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FIG. 3: The goodness of approximation F(r), denned in Eq. (49), for both the over- and under-density. The expansion up 
to 2nd order in gradients it plotted in solid red, up to 4th order in gradients in dotted blue. The units are irrelevant. What 
matters is the improvements of the 4th order over the 2nd order solution. See the text for an explanation of the shape of these 
curves. 



we show the solution at first order in linear perturbation theory, which shows no sign of collapse. The fourth order 
gradient expansion follows the exact solution closer than the second order expansion, but both anyway show the 
collapse to a central singularity. 

In order to quantify the improvement of the expansion up to fourth order in gradients over the expansion up to 
second order in gradients, we calculate the square of the difference between the approximation and the exact solution: 

f(r, t) = [X gradicnt (r, t) - A oxact (r, t)] , 

F(r)= f° f(r,t)dt. (49) 

Jo 

One could also integrate f(r,t) over time, or both time and radius, but we focus on the time integrated goodness of 
fit at different radii. This way we can see in Figure 3 that the expansion up to fourth order in gradients everywhere 
outperforms the expansion up to second order reaching a time integrated error yF of less than 10% for all radii. 
In both configurations there is a point at r ~ 0.3L where the first gradient term Rkij — ARij = which signifies 
the transition between the central region and the compensating shells. At this point the leading order correction is 
actually the 4-gradient term and thus offers the smallest improvement in approximation; the point behaves closest to 
a "separate universe" . A better approximation at this point would presumably require the inclusion of higher order 
terms. The overall approximation is better in the central region where the perturbation profile is relatively flat. In 
the compensating shells, r > 0.3i? the value of y/~F is higher but the inclusion of the 4-gradient term offers significant 
improvement. Eventually the approximation converges the the flat FLRW metric as r — > L. The reason why F(r) 
does not go to zero for r > L, is simply that we are comparing a numerical integration to an exact solution, and the 
difference between the solutions agrees within the desired numerical accuracy in the integration process, which is an 
arbitrary choice. 



IV. DISCUSSION 



The gradient expansion is an approximation scheme of a different nature from standard cosmological perturbation 
theory. Instead of assuming perturbations of small amplitude around a background, it builds a solution from an 
inhomogeneous initial seed metric kij by including terms with an increasing number of gradients of kij through 
combinations of the 3-Ricci tensor R\j\kij\ and its derivatives: %j ~ J2 ^n(t)(R r ij + $3 V 2r .R , ' -r ), where a hat 

n r 

indicates that quantities are computed from kij. The time evolution is then encapsulated in the functions F n (t). The 
approximation is expected to break down when F\(t)R ~ 1 but this roughly coincides with the time of collapse of 



\s P \ 



> 1 



over-dense regions. The ability of describing collapse and consequently the evolution in the regime when 
clearly goes beyond the capability of conventional perturbation theory at any order. The gradient series also contains 
the Zel'dovich approximation in its lowest order while higher orders provide relativistic corrections which are in fact 
required for greater accuracy. 

In this paper the gradient expansion was formulated for the case of a ACDM universe and offers a way to approximate 
structure formation in the non-linear regime. Furthermore, in order to gauge the accuracy of the approximation, we 
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have used it on a spherically symmetric distribution of matter with a non-trivial radial profile for which exact solutions 
are known in terms of the LTB metric. We found the approximation to work rather well. The first non-trivial terms 
containing 2 gradients describe the situation qualitatively well by showing collapse to a singularity in a manner 
comparable to the Zel'dovich approximation. The inclusion of 4 gradients provides quantitative improvement by 
approaching the exact solution even better and it is expected that the inclusion of higher order terms would provide 
even more accuracy, an expectation that would be interesting to confirm in future work. The expansion could also 
be checked for less symmetric situations, perhaps by comparing with other exact solutions or numerical simulations. 
Even if at the end of the day one has to resort to simulations for the most accurate description of the non-linear 
dynamics, it would pay to have analytic approximations that capture the essential features of non-linear evolution 
in a relatively simple way and with some level of realism as they would be applicable to a variety of problems. For 
this reason gradient-expansion techniques might be be useful additions to the toolbox of cosmologists as they strive 
towards a more accurate description of the universe. 
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